# generatedvh-test-singleslice.py

# Test code to generate a dDVH from a DICOM-RT file(s).
# Roy Keyes (roy.coding@gmail)
# 4 Mar 2010

import dicom
import numpy as np
import pylab as pl
import matplotlib.nxutils as nx


def gendvhexmpl(cont,ctslice):
    '''Example of dDVH generatation from RTSS and RTDOSE files.'''

    ss = dicom.read_file('rtss.dcm')
    rtd = dicom.read_file('rtdose.dcm')

    # Grab a set of contour data
    clist=ss.ROIContours[3].Contours[cont].ContourData

    # Convert contour to set of triples
    j = 0
    jmax = len(clist)
    ca=[]
    while j < jmax:
        ca += [clist[j:j+3]]
        j += 3

    ca = np.array(ca)

    # Dose data from pixel array
    dose = rtd.pixel_array[ctslice,:,:]

    # Make coordinate array
    xspace,yspace = rtd.PixelSpacing # mm spacing bewteen points
    ashape = dose.shape

    # Image position offset
    pos = rtd.ImagePositionPatient # mm

    ###NOTE####
    # Coordinate order is not uniform between dose, contour, offset, etc arrays
    # Therefore he we use dose x for contour y and vice versa...
    ###########

    xoffset = pos[0]
    yoffset = pos[1]

    coord = []
    j,k = 0,0
    while j < ashape[1]:
        k = 0
        while k < ashape[0]:
            coord += [[j*xspace + xoffset,k*yspace + yoffset]]
            k += 1
        j += 1

    coord = np.array(coord)
    print 'coord.shape:',coord.shape
    print 'coord.shape:',coord.shape
    print 'ca.shape:',ca.shape
    print 'coord[:,0].max():',coord[:,0].max()
    print 'coord[:,0].min():',coord[:,0].min()
    print 'coord[0] span:',abs(coord[:,0].max() - coord[:,0].min())
    print 'coord[:,1].max():',coord[:,1].max()
    print 'coord[:,1].min():',coord[:,1].min()
    print 'coord[1] span:',abs(coord[:,1].max() - coord[:,1].min())

    # Exclude z (transverse) dimension from contour array
    inside = nx.points_inside_poly(coord, ca[:,0:2])
    print 'inside.shape:',inside.shape

    # Unravel dose array to match coord array
    # Dose must first be transposed to reorder indices to match
    #   coordinate array index order
    dose = dose.T.reshape(ashape[0]*ashape[1])

    def filt(a,b):
        if a:
            return b
        else:
            return []

    # Filter out points outside of ROI
    roipoints = filter(lambda x:x,map(filt,inside,dose))
    #print 'roipoints:', roipoints

    # DEBUGGING: Create list of ROI filtered coords
    roicoord = filter(lambda x:len(x)>0,map(filt,inside,coord))

    return np.array(roipoints), np.array(roicoord), np.array(ca[:,0:2]),coord

# set contour and slice number
rp0,rc0,rpoly0,c0 = gendvhexmpl(0,0)
#rp1,rc1,rpoly1,c1 = gendvhexmpl(1,0)
#rp2,rc2,rpoly2,c2 = gendvhexmpl(2,0)
#rp3,rc3,rpoly3,c3 = gendvhexmpl(23,7)

print 'len(rp0):',len(rp0)
#print 'len(rp1):',len(rp1)
#print 'len(rp2):',len(rp2)

#print 'len(roipoints):',len(roipoints)

#n, bins, patches = pl.hist(roipoints, 50, normed=1, facecolor='green', alpha=0.75)

#roicoord = np.array(roicoord)
#print rc1

pl.plot(c0[:,0],-c0[:,1],'o')
pl.plot(rc0[:,0],-rc0[:,1],'x')
#pl.plot(rc1[:,0],-rc1[:,1],'x')
#pl.plot(rc2[:,0],-rc2[:,1],'x')
#pl.plot(rc3[:,0],-rc3[:,1],'x')

pl.plot(rpoly0[:,0],-rpoly0[:,1],'x')
#pl.plot(rpoly1[:,0],-rpoly1[:,1],'x')
#pl.plot(rpoly2[:,0],-rpoly2[:,1],'x')
#pl.plot(rpoly3[:,0],-rpoly3[:,1],'x')

pl.xlim(-229,254)
pl.ylim(100,420)

pl.show()


